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Abstract 


Zero-temperature density functional theory (DFT), density functional theory molec- 
ular dynamics (DFT-MD), and classical molecular dynamics using polarizable force 
fields (PFF-MD) are employed to evaluate the influence of Li + on the structure, trans- 
port, and electrochemical stability of three potential ionic liquid electrolytes: N— 
methyl- IV-butylpyrrolidinium bis(trifhioromethanesulfonyl)imide ([pyrl4] [TFSI] ) , N— 
methyl- IV-propylpyrrolidinium bis(fluorosulfonyl)imide ( [pyr 13] [FSI] ) , and l-ethyl-3— 
methylimidazolium boron tetrafluoride ([EMIM][BF 4 ]). We characterize the Li + sol- 
vation shell through zero-temperature DFT simulations of [Li(Anion) n ]( n_1 ) _ clusters, 
DFT-MD simulations of isolated Li + in small ionic liquid systems, and PFF-MD simu- 
lations with high Li-doping levels in large ionic liquid systems. At low levels of Li-salt 
doping, highly stable solvation shells having 2-3 anions are seen in both [pyrl4] [TFSI] 
and [pyrl3] [FSI] , while solvation shells with 4 anions dominate in [EMIM][BF 4 ]. At 
higher levels of doping, we find the formation of complex Li-network structures that 
increase the frequency of 4 anion-coordinated solvation shells. A comparison of compu- 
tational and experimental Raman spectra for a wide range of [Li (Anion) n ]( n_1 ' clusters 
shows that our proposed structures are consistent with experiment. We estimate the 
ion diffusion coefficients and quantify both size and simulation time effects. We find 
estimates of Li + diffusion are a reasonable order of magnitude and can be corrected 
for simulation time effects. Simulation size, on the other hand, is also important, with 
diffusion coefficients from long PFF-MD simulations of small cells having 20-40% error 
compared to large-cell values. Finally, we compute the electrochemical window using 
differences in electronic energy levels of both isolated cation/anion pairs and small ionic 
liquid systems with Li-salt doping. The single pair and liquid-phase systems provide 
similar estimates of electrochemical window, while Li-doping in the liquid-phase sys- 
tems results in electrochemical windows little changed from the neat systems. Pure 
and hybrid functionals systematically provide an upper and lower bound, respectively, 
to the experimental electrochemical window for the systems studied here. 
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Introduction 


Room temperature ionic liquids have been proposed recently as electrolytes for both 
conventional Li-ion batteries 1-3 as well as advanced high energy-density rechargeable bat- 
teries. 4 h In conventional Li-ion batteries, ionic liquids with organic additives have been 
shown to be safe and stable electrolytes with comparable performance to inherently more 
volatile pure organic electrolytes. 1-3 Other ionic liquid variants, particularly those containing 
imide-based anions, have been shown to form a passivating, Li + -permeable surface layers on 
the high energy-density electrode material, Li-metal. 4-6 Enabling the use of Li-metal would 
additionally aid the viability of advanced battery chemistries, such as lithium-sulfur and 
lithium-oxygen. ' ,8 Thus, it is important to understand the influence of Li + on the structure 
and properties of these potential electrolytes. In that direction, an abundance of recent ex- 
perimental and theoretical works have focused on the Li + solvation structure, 9-13 transport 
properties, 14-19 and electrochemical performance. 20-25 

In ionic liquids, the Li + solvation shell and conductivity are intimately related, with 
larger solvation shells generally having both a longer residence time and decreased ionic 
conductivity. 18,26 When the solvating anions are relatively complex, the resulting solvation 
structure may not be intuitive, such as with the [TFSI] and [FSI] anions (both of which are 
leading candidate anions for battery electrolytes. 4,6 ). Along these lines, a number of joint 
experimental-computational studies have been performed to elucidate the solvation shell in 
ionic liquids having the [TFSI] 11,13 and [FSI] 12 anions. Interestingly, spectroscopic mea- 
surements indicate that Raman-active modes related to the S-N and S-0 stretching modes of 
these anions are perturbed to higher frequency in the presence of Li + . Using the computed 
Raman activity of [Li(Anion) n ]( n-1 )“ clusters, one may identify which solvation structures 
correlate best with the perturbed frequencies and thereby determine the coordination num- 
ber, n. For [TFSI], the DFT Raman activities have been determined for n = 2 and n = 4 
clusters, with the n = 2 clusters 9-11 providing the best agreement with experimental Raman 
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spectra. For [FSI], on the other hand, higher coordinations seem to be preferred, 12 with the 
computed Raman activities of n = 3 clusters closely agreeing with experiment. Complemen- 
tary higher energy X-ray diffraction studies 27 on these ionic liquids suggest tetrahedrally 
coordinated Li + ions, and further lend credence to Li + being doubly and triply bound by 
[TFSI| and [FSI], respectively. 

As a property closely correlated with solvation structure, the diffusion coefficient of Li + 
is important when assessing the suitability of a given electrolyte. Experimentally, diffusion 
coefficients of Li + in ionic liquids are often obtained through pulsed-field-gradient spin-echo 
NMR. 2,15 However, experimentally providing a full measure of transport across a wide range 
of temperature and Li-salt mole fraction (xy) has only been attempted for a few ionic liq- 
uids, notably those having the [TFSI] and [BF 4 ] anions. 16,28 While such measurements are 
invaluable, MD simulations performed with polarizable force fields (PFF-MD) are an at- 
tractive alternative, having matured to the level where transport properties of neat ionic 
liquids 14 ’ 29 ~ 31 as well as those in the presence of Li-salts 17,26,31-33 can be determined within 
20-30 % of experiment. Accurately capturing the structure and energetics of the Li + sol- 
vation shell in a classical force held can be challenging. However, an alternative means of 
obtaining transport coefficients, though one still in its infancy, is through the application of 
density functional theory molecular dynamics (DFT-MD) simulations. Such simulations have 
recently been attempted for organic electrolytes, namely ethylene carbonate and dimethyl- 
carbonate. 34,35 Despite errors resulting from short simulation times and small sizes, results 
from such approaches are in reasonable agreement with experiment. For ionic liquids with 
Li-salt doping, however, the applicability of DFT-MD to transport is an open question, as 
the dynamics of ionic liquids are glassy at room temperature. This problem is further ex- 
acerbated at higher values of Li-salt mole fraction, x Li , where hundreds of picoseconds to 
hundreds of nanoseconds of simulation time can be required to obtain reasonable transport 
estimates. 26 

In addition to the molecular- level structure and transport properties, the electrochemical 
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window sets the limits of voltage bias reasonably accessible to a battery and is therefore 
important for assessing the performance of an electrolyte. The electrochemical window is 
generally experimentally estimated by characterizing current density as a function of volt- 
age . 36-38 Discontinuities in the current at high and low potential bias suggest electrolyte 
decomposition and set the anodic and cathodic limits, respectively, of the window. From a 
simulation standpoint, it is highly attractive to have a fundamental description of this quan- 
tity, particularly one that can be quickly determined as a screening tool , 39 and a number of 
theoretical approaches have recently been suggested to this end . 20-25 The most accurate, and 
complex, in this respect, determines the electrochemical window from the free energy change 
between reactants and potential decomposition products. To do this, a thermodynamic cycle 
is employed to connect the electrolyte free energy of decomposition in the gas phase to the 
fully solvated reactants in the liquid phase . 21 This procedure provides good agreement with 
experimental electrochemical windows of various ionic liquids, which include those having 
imidazolium and pyrrolidinium cations and imide and tetrafluoroborate anions . 24 Alterna- 
tively, reasonable experimental agreement has also been obtained by estimating stability 
from the electron affinity and ionization energy of individual molecules. This has been per- 
formed both in isolation and with solvent effects , 20-22 and does not require the determination 
of complex decomposition products. A final means of estimating electrochemical stability is 
from the difference between the highest unoccupied and lowest occupied unperturbed elec- 
tronic states, which has shown reasonable agreement with experiment and is computationally 
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inexpensive. ’ 

The present work is an attempt to apply in a broad and cooperative manner zero- 
temperature DFT, DFT-MD, and classical PFF-MD simulations to evaluate properties that 
characterize electrolyte performance. These include the previously discussed Li + solvation 
shell structure, the diffusion coefficient of Li + , and the electrochemical stability of the elec- 
trolyte as a function of xlu To benchmark these methods, we have here, as with our previous 
work characterizing energetic and transport properties , 26,40 considered three of the leading 
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ionic liquid electrolyte candidates for conventional and advanced batteries: [pyrl4] [TFSI] , 
[pyrl 3] [FSI ] , and [EMIMHBF4]. We begin our study with an investigation of the Li + sol- 
vation shell. Zero-temperature DFT is employed to analyze the structure and energetics of 
various [Li (Anion) n ]( n_1 D clusters. We then determine the Raman activities of the cluster 
using zero-temperature DFT, with an experimental comparison suggesting which clusters 
most closely correspond to liquid solvation structures. The stability of the most promising 
solvation shells are then explicitly tested using 20 ps DFT-MD simulations of a single Li + 
in ionic liquid systems having 8, 10, 12, and 16 ion pairs. In conjunction with the DFT-MD 
studies, room-temperature PFF-MD simulations of ionic liquids having 0.05 < x Tji < 0.33 are 
performed to quantify Li-networking effects on the solvation shells and validate experimen- 
tally derived anion coordination numbers, n. We then turn to an investigation of transport 
by extending select DFT-MD simulations to 100 ps to determine the ion diffusion coefficients, 
with particular focus on Li + . Furthermore, we benchmark size and simulation length effects 
inherent to the DFT-MD estimates of ion diffusion through the use of PFF-MD simulations 
of small ionic liquid systems. Finally, we cooperatively employ zero-temperature DFT and 
PFF-MD to evaluate the electrochemical window in both neat ionic liquid electrolytes and 
those having Li-salt. In this respect, zero-temperature DFT is employed to evaluate both 
the energetic difference between the lowest unoccupied and highest occupied molecular or- 
bitals as well as the difference between ionization potential and electron affinity of isolated 
cation/anion pairs. In the liquid phase, PFF-MD is used to generate room-temperature 
structures of neat and Li-doped ionic liquids having 24 ion pairs, the electronic structure of 
which was determined to be free of size effects, that are then combined with DFT compu- 
tations to provide averaged values of the energetic difference between the conduction band 
minimum and valence band maximum. We perform these computations with a mixture of 
both pure and hybrid functionals, which for our 24 ion pair simulations represent amongst 
the largest systems treated with these methods. 
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Methods 


Clusters: DFT Computation 

As a fundamental means of understanding the Li + solvation structure in the ionic liquids 
of interest here, we perform various DFT computations of isolated [Li (Anion) n ]( n_1 ' clusters, 
where 2 < n < 4. For each individual n-state, we have performed a detailed study of cluster 
energetics with respect to configuration, systematically probing structures with different 
combinations of single-fold, monodentate bonds ( 77 1 ) and two-fold, bidentate (■ r f) ligands. 
For each cluster, the harmonic vibrational frequencies are computed and confirmed to be 
positive to ensure that the clusters represent stable minima. Additionally, for each harmonic 
frequency we compute the associated IR intensity and Raman activity. Because many of the 
clusters are too large to treat with high-level quantum chemical methods and large basis sets, 
we have performed these calculations with the B3LYP functional 41,42 in combination with 
the 6-31+G** basis set of Pople and coworkers . 43 In a previous work, we showed that this 
combination of theory and basis set provides reasonable agreement with second order Mpller- 
Plesset perturbation theory in the complete basis set limit for both Li + /anion interactions 
and for n = 2 clusters . 40 We have performed these computations with Gaussian 09. 44 

Additionally, we have investigated the electrochemical stability of our ionic liquid elec- 
trolytes through a DFT analysis of isolated cation/anion pairs. The focus of our compu- 
tations were two-fold: understanding the influence of functional on estimated stability and 
a comparative assessment of multiple measures of electrochemical stability. With regard 
to the former, we have employed a wide variety of both pure exchange/correlation func- 
tionals, which include the exchange functional of Becke 45 with the correlation functional of 
Perdew and Wang (BPW91), the exchange and correlation functionals of Perdew and Wang 
(here termed PW91), the exchange and correlation functional of Perdew, Burke, and Ernz- 
erhof (here termed simply PBE ), 46,47 and various hybrid functionals including the screened 
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Coulomb functional of Heyd, Scuseria, and Ernzerhof (HSE06), 48,49 B3LYP, and M06 from 
Truhlar and coworkers. 50 With regard to the latter, we have compared both the difference 
in the lowest unoccupied molecular orbital (LUMO) and highest occupied molecular orbital 
(HOMO) energies (AEhl) as well as the sum of ionization potential (IP) and electron affinity 
(EA, which is here defined as the energy released upon addition of an electron), AEip, to 
experimentally evaluated electrochemical windows as an assessment of their relationship to 
electrochemical stability. As with the [Li(Anion) n ] _n+1 clusters, we have chosen to perform 
these computations with Gaussian 09. 44 

Liquid Simulations: Ab Initio Molecular Dynamics 

To extend the cluster computations to include the effects of temperature and explicit 
solvation, we have performed a set of 100 ps DFT-MD simulations on 7[pyrl4]|TFSI] + 
Li[TFSI], 9[pyrl3] [ESI] + Li[FSI], and 11[EMIM][BF 4 ] + Li[BF 4 ] systems at T = 363 K, 
at which temperature the ion mobility is such that multiple Li + solvation variations can 
be sampled within the stated 100 ps timeframe. These simulations are performed with the 
Vienna Ab Initio Simulation Package (VASP) 51-54 using the frozen core all electron projector 
augmented wave (PAW) method 55,56 in the generalized gradient approximation of Perdew, 
Burke, and Ernzerhof. All DFT-MD simulations are Gamma-point computations with an 
energy cutoff of 400 eV, an electronic energy convergence criteria of lxlO -4 eV, and a time 
step of 0.5 fs. 

Along a similar vein to the previously described electrochemical stability study on isolated 
cation/anion pairs, we have additionally carried out computations to determine the difference 
between the conduction band minimum (CBM) and the valence band maximum (VBM), 
given as AEyc, of liquid phase ionic liquids using DFT and a plane- wave basis set. In this 
case, we initially perform classical molecular dynamics simulations, which will be described 
in detail later, at T = 298 K on systems with up to 24 ionic liquid pairs to obtain trajectory 



information over a 6 ns simulation length. Using the trajectory information so obtained, we 
then choose 10 representative configurations throughout the simulation and, by performing 
single point DFT computations on each configuration, estimate a temperature average of 
AEyc, which is akin to the value of AEhl for single pairs. Using the largest 24 ionic 
liquid pair system, we then replace cations with Li + to produce systems with Li-doping 
levels up to xl; = 0.33. The same coupled classical molecular dynamics-density function 
theory approach for estimating AEyc in neat ionic liquids is again performed with each of 
these systems to produce a temperature and explicit solvent estimate of the electrochemical 
window. As with AEhl, AEyc from liquid phase computation is computed employing both 
the pure functionals PBE and PW91 as well as the hybrid functionals HSE06 and B3LYP. We 
note that fully 3-D periodic computations with hybrid functionals having many thousands 
of electrons are computationally expensive and represent the upper size limit tenable with 
these methods an our computational resources. 

Liquid Simulations: Polarizable Molecular Dynamics 

Complementary to the first principle MD simulations, we have performed a suite of liquid- 
phase molecular dynamics simulations using the atomistic polarizable potential for liquids, 
electrolytes, and polymers (APPLE&P), which has been shown to offer unparalleled accu- 
racy in both the structure-energy relationship and transport of ionic liquids in the presence 
of Li-salts. 14 ’ 17,26 ’ 30 ’ 31 ’ 33 ’ 40 ’ 57 ’ 58 The polarizable force held MD (PFF-MD) simulations are 
performed with a modified version of the large-scale atomic/molecular massively parallel 
simulator (LAMMPS), 59 with the specific details of the polarizable force held implementa- 
tion and MD simulations for these liquids being found in a previous work. 26 While there are 
multiple purposes for the use of classical molecular dynamics in the simulation of ionic liq- 
uids, our primary purpose here is to use PFF-MD as a platform for understanding potential 
system size influences on the DFT-MD results. As such, throughout the present work we 
provide multiple comparisons of polarizable molecule dynamics simulations (PFF-MD) with 
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measures of structure and transport derived from DFT-MD. 


Results 

Li + Solvation Structure: Cluster Computations 

Energetics 


To initiate our study of Li + solvation by ionic liquid, we have probed the configuration 
space of [Li(Anion) n ]0 _1 ) _ clusters, with the number of anions, n, being varied from one 
to four, with DFT computations at the B3LYP/6-31G** level. As previous theoretical and 
experimental 9,10,12,27 ’ 60-64 studies have suggested the number of Li-solvating anions for the 
liquids of interest here to be 2 < n < 4, we herein limit our analysis to clusters with n < 4. In 
a previous work, we have found that in Li|BF 4 ] clusters with n = 2 or 3, if bonding dominates 
and produces the most energetically stable clusters at T — O. 40 Upon reaching clusters with 
n = 4, the energy of bringing an anion close enough to form a if bond is outweighed by 
the mutual anion-anion repulsion, which results in a stable 4 if cluster. With this in mind, 
we have extended our previous work to include n = 2, 3, and 4 Li + clusters with both the 
cis and trans conformers of the [TFSI] and [FSI] anions, the results of which, along with 
|BF 4 ], are summarized in Table 1, with the structures associated with the most energetically 
favorable clusters being shown in Figure 1. To simplify the following discussion, we note 
distinct trends in the conformational preference of a [TFSI] or [FSI] anion in a given cluster 
with Li-binding. It was found that clusters having trans- 1 TFSI] are the most energetically 
stable, while [FSI] is more complex, with rf anions taking the cis conformation and rf anions 
taking the trans conformation. Only the low energy clusters following these rules are shown 
in Table 1, and a full survey of the energetics of more unstable clusters that deviate from 
these rules is shown in the supplementary. 
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To elaborate upon the energetics of the clusters, we have decomposed the energy into the 
primarily repulsive interactions of the anions (Er), which is the energy of bringing the anions 
into their cluster configuration from infinite separation without Li + ; the primarily attractive 
interaction between the Li + and the anions (Ea), which is obtained from the difference in 
total energy and Er; and the cluster binding energy (Er, where Er = — (Er + Ea). In 
agreement with previous computations, 9 ’ 10,12 ’ 40 the n = 2 clusters assume rj 2 configurations, 
while the n = 4 clusters assume 77 1 configurations. For the n = 3 clusters, all anions evaluated 
exhibit higher values of Ea with increased r f binding. In fact, as can be seen from Table 1, 
the conversion of any r/ 1 anion to an r / 2 state decreases Ea by roughly 8-12 kcal/mol across 
all anions investigated. On the other hand, Er exhibits the converse relationship, with more 
tf binding leading to higher levels of repulsion, which is on the order of 8-11 kcal/mol for 
each z/ 1 anion converted to an rj 2 state. For n = 3 clusters with [TFSI] and [FSI] anions, 
these energetic contributions appear to provide a minimum Er in the r) 2 2r) 1 cluster, being 
4.3 kcal/mol and 2.7 kcal/mol lower in energy than the fully bidentate 3r] 2 state in [TFSI] 
and [FSI], respectively. A comparative analysis of Er between clusters with different values 
of n is difficult, as higher values of n necessarily lead to lower values of Er due to increased 
anion repulsion. Solvation effects in a charge neutral system must be considered to provide 
more clarity on the relative stability of the Li + solvation shell with respect to n. 

Raman Analysis 

To experimentally corroborate the structure and number of anions in the Li + solvation 
shell, we now turn to a study of Raman spectroscopy. For [TFSI] and [FSI] anions, Raman 
spectroscopy allows the identification of characteristic anion vibrations that fundamentally 
change when bound to Li + . By knowing the ratio of Li + to anions, the ratio of the magnitude 
of these modes can be directly related to Nj, the total number of anions solvating all Li + , 
as well as n, the characteristic number of anions in the solvation shell. With n at hand, 
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the details of the solvation structure can be predicted through DFT cluster computations 
akin to those presented in Table 1. At room-temperature and values of xl; < 0.2, such 
procedures have predicted Li + to be coordinated at room temperature by 1.86-2 [TFSI] 
anions in a 2 if configuration 9 11 and 2.9-3 [FSI] anions in a rj 2 2r) 1 configuration, which 
reaches an equilibrium with an n = 2, 2 r] 2 structure at higher temperature. 12 While our PFF- 
MD and DFT-MD simulations do not dispute these results, we additionally find that n = 3 
solvation structures are a present in [pyr 14] [TFSI] . To explore the feasibility of this structure 
and assess the influence of structural variation on Raman activities, we presently re-evaluate 
the available experimental data by performing a thorough comparison of the DFT-mediated 
Raman spectra of [Li(TFSI) n ]^ n_1 ^ _ and [Li(FSI) n ] ( n— c l us ters having 2 < n < 4. While 
the experimental hirelings have indeed also been aided by DFT Raman spectra, 9-12 we have 
here provided a more in-depth study on the influence of structure and value of n, especially 
for | TFSI] where n = 3 solvation structures have not been vibrationally characterized. 

In |TFSI|-based ionic liquids, Li-salt doping gives rise to a new Raman signature between 
745-750 cm -1 , which is a perturbation of the coupled CF3 bend, 5 S (CF3), and SN stretch, 
u s ( SN), present at 742 cm -1 in neat samples. The experimental Raman spectra of these 
modes in neat |BMIM][TFSI| and [BMIM][TFSI] with x Li = 0.33 Li|TFSI], as taken from 
the work of Lassegues and coworkers, 9 is displayed in Figure 3d. As can be seen, when Li-salt 
doping levels increase, the higher frequency signature from [TFSI] bound to Li + increases, 
while the lower frequency signature, which is assumed to be entirely [TFSI] not bound to 
Li + , decreases. Along with the reproduction of the experimental data in Figure 3, we have 
included our computed Raman activities for ]Li(TFSI) n ]( n_1 ' clusters having n = 2, 3, and 
4. The frequencies and activities shown were computed with B3LYP/6-31+G** for the 2 rf 
and r) 2 2r) 1 structures that we find using DFT-MD as well as the Ar/ 1 structure that, while not 
observed through DFT-MD, has been previously proposed as a potential solvation structure. 
To provide a more thorough estimate of Raman activities, we have additionally computed 
these structures in both their cis and trans conformations as well as with multiple local 
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geometries, which in total lead to 2, 6, and 10 individual clusters computations for the 2 if, 
r] 2 2i] 1 , and ir/ 1 clusters, respectively. As can be seen in Figure 3, n = 2 and n = 3 cluster 
have Raman active modes in the 745-755 cm -1 range, while n = 4 clusters appear to have 
modes shifted to lower frequencies in the 735-745 cm -1 range. While n = 2 has a strong 
activity only in the high-frequency region, the n = 3 clusters have activities divided between 
the low and high frequency regions, which is problematic when assigning all Li-bound anions 
to the high-frequency region. From further analysis of our computed vibrational modes, it 
appears that primarily [TFSI] having r/ 2 binding appear in the high-frequency region, while 
77 1 is indistinguishable from the neat frequency. Therefore, given that experimental analysis 
associates 1.8-2 anions per Li + as coming from the 745-750 cm -1 signature, a potential 
scenario is a coexistence of n = 2 and n = 3 clusters. 

A similar Raman analysis can be performed for [FSI], in which both the i/ s (SN) at 731 
cm -1 and the z/ s (SO) at 1220 cm -1 exhibit a strong dependence on Li-salt doping. As 
with |TFSl|-based ionic liquids, increasing the Li-salt doping in an [FSI|-based liquid results 
in a higher- frequency emission, which occurs at 744 cm -1 for the z/ s (SN) vibration and 
1230 for the z/ s (SO) vibration, associated with anions bound to Li + progressively overtaking 
the lower-frequency emission of the free [FSI]. Experimental Raman spectra displaying the 
z/ s (SN) and y s (SO) Raman-active modes, as taken from Fujii and coworkers, 12 are shown in 
Figures 4d and 5d, respectively, for [EMIM][FSI] in both the neat form and that having Xl; = 
0.225 Li[FSI] . The associated DFT-derived Raman activities in these regions are included in 
Figures 4 and 5 for n = 2, 3, and 4 clusters. Similarly to the case of [TFSI] clusters, we have 
performed computations on clusters with cis and trans conformers as well as with multiple 
local geometries, again leading to 2, 6, and 10 individual clusters computations for the 2 if, 
r) 2 2r) 1 , and 4 r/ 1 clusters, respectively. Concerning the z/ s (SN) vibrations, [FSI] associated 
with n = 4 clusters has Raman-active modes primarily appearing in the neat region, while 
[FSI] associated with n = 2 and 3 clusters has higher frequency Raman active modes. For 
the za.(SO) mode, we see a slightly different behavior, with n = 3 and n = 4 clusters having 
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activities arising in the high-frequency region and n = 2 clusters having activities appearing 
in the low-frequency region. In contrast to [TFSI], some t] 1 z/ s (SN) vibrations in [FSI] appear 
in the high-frequency regions, which leads to n = 3 clusters having frequencies spanning a 
large region between 730 and 770 cm -1 . These activities are, however, roughly centered 
around 744 cm -1 , which corresponds to the experimental Raman-active frequency of [FSI] 
bound to Li + . The u s (SO) mode of [FSI] appears to be more susceptible to Li-binding as most 
configurations of n = 3 and 4 clusters have strong activities in the high-frequency regions, 
with the exception of n = 2 clusters that, anomalously, appear at low-frequency. Taking this 
with the experimental, room-temperature coordination of 2.9-3, we must conclude that n = 3 
clusters are indeed dominant, and, further, at high temperatures, where coordination number 
via the z/ s (SO) mode is found to decrease to 2.5, the loss of activity in the high-frequency 
region can ascribed to the evolution of more n = 2 clusters. 

Li + Solvation Structure: Molecular Dynamics 

DFT-MD Stability of Solvation Shells 


A natural extension of the [Li(Anion) n ]( n- h“ cluster study is to introduce both tem- 
perature effects and solvation effects, which we accomplish here by performing liquid-phase 
DFT-MD simulations of our electrolytes in the presence of Li + . To do this, for all three ionic 
liquids we have employed PFF-MD simulations to generate small, liquid-phase simulation 
cells having 8, 10, 12, and 16 pairs, with one cation in each system being replaced by Li + . 
The structures are initially equilibrated using PFF-MD at 363 K for 10 ns, and multiple 
configurations of Li + having n = 2, 3, and 4 solvation shells are chosen as starting states 
for DFT-MD simulation. We perform 20 ps DFT-MD simulations with the system sizes 
and configurations so described, and we find that for the elevated T = 363 K simulation 
temperature is sufficiently long to sample multiple solvation shell configurations. 

For [pyrl4] [TFSI] and [pyrl3] [FSI] , we fold n = 2 and 3 solvation structures, as shown 
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in Figure 1 a and b, to be stable over the 20 ps DFT-MD simulations, with the n = 4 
solvation shells rapidly converting to an n = 3 state. Moreover, we find that the preferred n 
= 3 solvation shell corresponds to the most stable r) 2 2r) 1 cluster as evaluated from Table 1. 
On the other hand, the solvation shell of Li + in |EMIM]|BF 4 | is primarily 4r/ 1 , as shown in 
Figure 1 f, with the solvation shell transitionally assuming the r) 2 2r) 1 configuration, which 
agrees with the energetic trends noted in Table 1. From the present set of simulations, we 
have not observed strong size-related effects on the solvation shell, with our smallest systems 
producing the same combination of solvation shells as our largest. These results are also 
in generally good agreement with the experimentally supported n = 2 [TFSI] 10,63 and n = 
3 [FSI ] . 12 However, the n = 3, |TFSI]-coordinated solvation structure has not heretofore 
been described by experiment; however, we suggest that such a structure is viable and is 
additionally asserted its through the previous comparison of computational and experimental 
Raman spectra. 

Expanding on our DFT-MD study, we have benchmarked the structural similarity of 
PFF-MD simulations to that of the DFT-MD simulations. To improve the statistics of our 
first principles molecular dynamics, we have performed an additional set of 100 ps DFT-MD 
simulations at 363 K beyond the initial 20 ps equilibration on the 8, 10, and 12 pair systems 
of [pyrl4] [TFSI] , [pyr 1 3] [FSI] , and [EMIM]|BF 4 ], respectively, which represent the largest 
systems that can be reasonably simulated with DFT-MD at such a time-scale. Throughout 
these long simulations, no new solvation structures were observed. In Figure 2, we com- 
pare the radial distribution functions of these systems as obtained with DFT-MD to those 
produced from PFF-MD simulations at the same size and temperature conditions. Across 
the board, the agreement between the radial distribution functions is excellent, especially 
when considering the smaller, 100 ps length of the DFT-MD simulations as compared to 
the 6 ns used for the PFF-MD simulations. In the [pyrl4] [TFSI] and [pyr 1 3] [FSI | cases, 
the initial g(r) peak centered about r = 3 A shows a double hump that is characteristic 
of the more distant r/ 1 anions and the closer rj 2 anions. For [BF 4 ], there is a much more 
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localized signature around 3 A, which encompasses both r) 1 anions and the rarer r] 2 anions. 
As signified by the agreement in both the shape and magnitude of the first g(r) peak, the 
solvation structures observed from the PFF-MD simulations are very similar to those seen in 
DFT-MD. However, these results also suggest that PFF-MD slightly under-predicts the net 
amount of rf binding compared with the DFT-MD simulations, especially for [TFSI] and 
[FSI] anions. 

PFF-MD Network Analysis 


While the previous analysis has been performed for the low Li-salt doping regime, we 
now discuss the case of high Li-salt doping and the potential influence of Li-networks on 
the experimental estimation of n. As Li-salt doping increases, the experimentally derived 
solvation number exhibits a marked decrease, with |TFSI|-based liquids decreasing to roughly 
n = 1.5 at xu > 0.15. 9 Such a decrease cannot be attributed to the solvation numbers of 
individual Li + , where coordination by one anion is not sensible, and must be related to the 
effect of Li-networks, where multiple Li + are bridged by mutual anions. Because current 
experimental practices can only determine total number of anions bound to all Li + , the 
formation of networks, where a given anion can belong the solvation shell of multiple Li- 
ions, leads to an underestimation of n. To quantify this effect we have performed 216 ionic 
liquid pair room temperature PFF-MD simulations having Xl; values of 0.05, 0.15, and 
0.33 and have computed the average size of Li. ..Li network, Mu... Li, the average number of 
anions in each individual Li + solvation shell, (n ) , and the total number of unique anions 
solvating all Li + , (Nj), divided by the total number of Li + , Nll The value of (N“)/Nl; is 
more aligned with that obtained from spectroscopic studies and is independent of the local 
solvation number of each anion, which simulation suggests to alway be 2 or 3. As shown 
in Table 3, the value of NLi...Li increases from 1.1-1. 2 at xl; = 0.05 to 1.5-1. 7 at xli = 0.33, 
indicating a shift from isolated Li-ions at low Li-doping levels to Li-networks at moderate 
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to high Li-doping levels. Furthermore, we note an inverse relationship between (n) and 
(N7 )/Nli with increasing x^; from = 0.05 to xl; = 0.33, we see a slight increase of 0.1- 
0.3 in (n ) , while (Nj)/NLi is decreased by one. The magnitude of Li-networking is notable 
in even the xl; = 0.15 system, where (Nj)/NLi decreases by a half. These results are in 
generally good agreement with experiment and suggest an under-prediction of n at high xl;; 
however, clarifying the exact magnitude of this effect would require a more quantitative look 
at the structure and vibrational frequencies of Li-networks, which is out of the scope of the 
present work. 

IR Spectra 


Whereas the previous discussion of Raman spectra focused on how Li + alters the native 
vibrations of anions in its solvation shell, we now turn to a brief discussion of the spectroscopic 
signature that arises from the vibration of Li + itself. Such vibration are most apparent at 
lower frequencies (< 500 cm -1 ) from IR spectroscopy, where IR active modes arise due to the 
vibration of Li + within its solvation shell, or the so-called “rattling" motions. In this respect, 
an interesting approach to characterizing the Li + IR signature is to leverage our long DFT- 
MD simulations at 363 K to compute the IR spectra. IR spectra computations from MD 
simulations are not commonly undertaken because the DFT-MD simulation time required 
to accurately resolve a frequency can be large for low-frequency vibrations and classical MD 
simulations are generally not sufficiently accurate to produce viable spectra. However, the Li 
“rattling" motions produce a broad IR signature that can be qualitatively characterized over 
the 100 ps time scale of our DFT-MD simulations. To do this, we first characterize multiple 
configurations from our DFT-MD simulations with Bader analysis to obtain an approximate, 
average charge for each atomic species. The net charge so obtained can then be coupled with 
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velocity data from the DFT-MD trajectory and used in an autocorrelation function, J, 

J(t) = <M(i) ■ M(0)) , (1) 

where M(t) = Y2l=o 'LwW- By simply taking the Fourier transform of J(t), one can then 
obtain the temperature-dependent IR spectrum of a given system, J(oj). 

The results of our DFT-MD procedure for our three ionic liquids are shown in Figure 6 and 
compared with the IR signature of the Li + vibrations obtained using DFT and the harmonic 
frequencies of all the clusters presented in Table 1. For the most part, the spectra computed 
from DFT-MD simulation are in qualitatively good agreement with the IR spectra produced 
from cluster simulations. Interestingly, independent of ionic liquid the Li + IR active modes 
are broadly spread about 300 cm -1 , which is in good agreement with experimental measures 
having a broad absorption around 374 cm -1 . 9 Interestingly, this IR signatures of the modes 
appear in the same frequency region for all three ionic liquids, which suggests the curvature 
of the Li-binding is similar in spite of the different preferred solvation shells and anions. 
We also note that the temperature-induced broadening from our DFT-MD simulations is 
well-approximated by the use of the many clusters present in Table 1, which, along with the 
previously described Raman study of cluster coordination number, suggests that multiple 
configurations of solvation complexes are required to obtain a reasonable measure of IR 
absorption or Raman activity. 

Transport Properties 

We now turn to an investigation of Li + transport from our 100 ps DFT-MD simulations. 
With regard to ion transport, PFF-MD using APPLE&P has proven an invaluable tool, with 
computed transport coefficients for ionic systems, which are traditionally inaccurate from 
classical MD simulations, being in close agreement with experimental measurements. 17 > 26 > 31 > 33 
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As an alternative to PFF-MD, one may employ DFT-MD simulations 34,35 in the liquid phase 
and determine the transport coefficients from the resulting trajectory data. However, due to 
the extraordinarily large computational cost of DFT-MD, there will ultimately be limitations 
on both the simulation size as well as the simulation duration, both of which strongly play 
into the statistical accuracy of the diffusion coefficient. To this end, we here benchmark the 
potential system size and simulation time effects present in our 100 ps DFT-MD simulations 
through the analysis of trends in equivalently sized PFF-MD simulations. As a means of 
understanding the simulation time influence on our DFT-MD simulations, we compare the 
apparent diffusion coefficients (given as D app (f), where t is time) of cations, anions, and Li + 
obtained from both our DFT and PFF simulations. The D app (f) values are a convenient 
means of understanding how error in the diffusion coefficient changes with simulation time, 
with the thermodynamic diffusion coefficient (D) being equal to the long time limit of the 
apparent diffusion coefficient, D“ = lim^oo D"’ app (f), where a is a particular ionic species. 
Both coefficients are, of course, related to the mean square displacement of the ionic species, 


N“ N” 


D“ = lim — / / 

t-¥ OO 6 t 

i = i i = i 


r a,t 


(i - S,) - r“-‘(0)] 2 > /(N“N°), 


( 2 ) 


where r Q,! is the atomic position of one of the N" atoms in molecular species a and tj is the 
j th time origin, with the diffusion being an average over N° time origins. 

In Figure 7, we present measures of D app (f) as a function of time for the various species 
in the ionic liquids of interest at 363 K. For the DFT-MD measures, given in Figure 7a- 
c, the total simulation times are 100 ps total simulation time, N° is taken as 20, and the 
apparent diffusion coefficient is shown up to t — 50 ps, which results in generally well- 
behaved D app (f) profiles. We note that diffusion within the relatively short 50 ps time- 
scale appears to still be within the ballistic regime; that is, the apparent diffusion is still 
correlated to its previous states and provides an overestimate of thermodynamic diffusion, 
with the thermodynamic limit being expected at large values of t where D app (f) is constant. 
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Despite this overestimate, the values of D app (t = 50ps) are the correct order of magnitude 
when compared to experimental 15-17 and previous computational 26 measures of diffusion and 
further follow the trend where the diffusion of cations (D + ) > anions (O') > Li + (D Li ). The 
values of D app (f) from PFF-MD using simulation cells equivalent to those of DFT-MD, as 
shown in Figure 7d-f, are well-converged and smooth, with the total simulation time being 
6 ns, N° being 300, and the apparent diffusion being shown up to t — 500 ps. These, too, 
follow the trend of D + > D~ > D Ll and agree in magnitude with the DFT-MD results. The 
much longer time-scales in PFF-MD simulations reveal that the ballistic regime is long-lived 
and on the order of hundreds of ps. 

For a more quantitative comparison of the diffusion estimates, we provide in Table 2 
measures of D app (f) as obtained from small DFT-MD and PFF-MD simulation cells at t — 
50 ps and for small and large PFF-MD simulation cells in the thermodynamic limit, given as 
t — > oo. First, comparing small PFF-MD and DFT-MD cells, we note that with only a few 
exceptions the DFT-MD values of D app (f = 50ps) are somewhat larger than those produced 
by PFF-MD; however, considering the error in the DFT-MD estimates, which is ~50% for 
the cations and anions and ~30% for the Li + ions, most of the PFF-MD results are within 
the error of the DFT-MD. To elaborate upon the diffusion of Li + in particular, we note the 
ratio of D Ll,app (t — > oo) to D Ll,app (t = 50ps) for the small PFF-MD systems is 0.37, 0.48, and 
0.77 for [pyrl4]|TFSI], [pyrl3] [FSI] , and [EMIM]|BF 4 ], respectively. Assuming the DFT-MD 
simulations would exhibit a similar decrease as t — > oo leads to Li + diffusion values of 0.71, 
0.78, and 0.8 x 10 10 m 2 /s, for [pyrl4] [TFSI] , [pyrl3] [FSI] , and [EMIM[BF 4 ], respectively. 
This is in good agreement with available experimental measurements (extrapolated to 363 
K) for [pyrl4] [TFSI] , 15,1 1 which range from 0.85-1.91 x 10 10 m 2 /s, and for [EMIM]|BF 4 ], 16 
which is estimated to be 0.8 x 10 10 m 2 /s, while the lack of experimental data for [pyrl3] [FSI] 
makes a meaningful comparison difficult. 

To quantify potential size effects related to the small-cell DFT-MD and PFF-MD results, 
we compare our small PFF-MD results to previously reported diffusion coefficients 26 as 
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computed from PFF-MD simulations having 144 ion pair for |pyrl4] [TFSI] and 216 ion pairs 
for both |pyrl3] [FSI] and [EMIM]|BF 4 ], To maintain a uniform value of xls across simulation 
size, we interpolate our large cell diffusion coefficient to the appropriate small cell value of 
XL, from values taken from independent simulations at 0.05, 0.10, and 0.15. The small-cell 
PFF-MD results are within 20-40% agreement with the large-cell results. For the time scales 
accessible to DFT-MD (~100 ps), the results suggest that the influence of simulation time, 
which can result in diffusion coefficients several times too large, is more important than 
simulation size, which results in 20-40% variance of diffusion coefficient. It is not clear, 
however, that the secondary influence of system size is a universality for ail system, with 
even previous MD simulations of other ionic liquids suggesting significant size effect up to 
many hundreds of pairs. 65 

Electrochemical Stability 


Single Ion Pairs 

As it currently stands, the majority of computational approaches for evaluating electro- 
chemical windows are based on electronic energy levels of isolated or solvated electrolyte 
molecules. 21-23,25 While more complex treatments exist, 20,24 we here follow and evaluate 
the validity of this convention on the basis of DFT-derived electronic energy levels. 21 23,25 
For isolated electrolyte molecules, two natural measure of electrochemical stability are the 
energy difference between the LUMO and HOMO, AEhl, and the difference between the 
electron affinity (EA) and ionization potential (IP), given here as AEiE,d for fully geometry 
relaxed, diabatic measures and AEie, v for instantaneous, vertical measures. An extension of 
these procedures can be made to the liquid phase, where the difference between HOMO and 
LUMO energies in isolated electrolyte molecules can be equated to the difference between the 
conduction band minimum (CBM) and valence band maximum (VBM) energies in periodic 
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DFT simulations. 


Table 4 displays estimates of electrochemical window on the basis of AEhl, AEie,<i, and 
AEie, v for single cation/anion pairs. For comparative purposes, we have presented results 
from a wide range of pure exchange/correlation functionals, including BPW91, PW91PW91, 
and PBEPBE, as well as a variety of hybrid functionals, including HSE06, B3LYP, and M06, 
which are all evaluated with the 6-31+G** basis set. Concerning AEhl, there is a clear 
divide between the pure and hybrid functionals, with the hybrid functionals yielding 1.5-2 
eV higher measures, in general agreement with the underestimation of the HOMO/LUMO 
gap in pure functionals. Similarly, for measures of AEiE,d and AEie, v hybrid functionals are 
0. 5-1.0 eV higher than the pure functionals. For almost all cases AEhl < AEiE,d < AEie, v , 
with AEhl being 3-4 eV lower than the measures from ionization potential and electron 
affinity; the exception to this is [pyrl3] [FSI ] as treated with pure functionals, where the 
[FSI] decomposed in the diabatic regime. Of interest is the fact that the value of EA is very 
small compared to IP for the liquids investigated, which suggests, from a vacuum reference, 
adding an electron to the cation is easier than removing one from the anion. 

Comparing the measures between different cation/anion pairs, the magnitude of AEhl 
between the three liquids are similar, with the variation between liquids for a given func- 
tional being < 0.6 eV. Similar trends can also be noted for measures of AEiE,d and AEie jV , 
though the variation between liquids is closer to 1 eV and values follow the general trend 
[EMIM||BF 4 ] > [pyrl3] [FSI] > [pyrl4] [TFSI] . The experimentally measured electrochemi- 
cal windows are 3.8-6 V, 36,37,66 5.3 V, 38 and 4.3 V 67,68 for [pyrl4][TFSI], [pyrl3] [FSI] , and 
[EMIM||BF 4 ], respectively. Generally speaking, the values as given suggest the IP+EA 
procedure leads to a gross overestimation of the electrochemical window, while the LUMO- 
HOMO procedure provides a lower and upper bound to the electrochemical window, with 
pure functional being the lower bound and the hybrid being the upper. The exception 
to this is [EMIM]|BF 4 ], where AEhl from pure functionals provides close agreement with 
experiment. 
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Liquid-Phase Simulations 


Following the implication from the single ion pair DFT computations that AEhl is a straight- 
forward means of bounding the electrochemical window, we now investigate the influence of 
temperature and solvation effects on these measures through coupled PFF-MD and DFT 
simulation. In this respect, we use PFF-MD as a platform for generating structures at room 
temperature and perform single point, DFT computations on select structures to generate 
an average measure of the difference between the CBM and VBM, given as AEyc, which is 
akin to the value of AEhl in isolated molecules. An alternative approach would be to employ 
the polarizable continuum medium approach on DFT-simulations of single ionic liquid pairs. 
However, it has recently been shown that the quadrupolar interactions present in the liquid 
solvent surrounding a target solute are difficult to approximate with continuum models, 20 
which, especially considering the ionic nature of our electrolyte, has led to our use of full 
liquid systems for the electrochemical window computations. As previously noted, 25 a seri- 
ous concern when performing liquid-phase simulations is the influence of system size on the 
resulting AEyc- To better understand potential size effects, we have performed PFF-MD 
simulations of liquid-phase systems having 12, 18 and 24 ionic liquid pairs, the simulation 
lengths of which are 60 ns at room temperature. From these simulations we have taken 10 
representative trajectories, which are 6 ns apart to ensure that the structures are signifi- 
cantly different, and performed electronic structure computations with DFT using the pure 
functionals PBEPBE and PW91PW91 and the hybrid functionals HSE06 and B3LYP. The 
results, shown in Table 5, indicate that for all liquids there is a significant structural influ- 
ence on the electronic energy levels, and thereby the estimate of electrochemical window, for 
systems numbering fewer than 18 ion pairs. 

Guided by our size effect study, we have made estimates of the electrochemical window 
of our electrolytes in both the neat form and with Li-salt doping. We have used our largest 
24 ion pair system and replaced cations with Li + to achieve systems having xli values of 
0.083, 0.166, 0.25, and 0.333. For each system, we have employed the previously described 
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procedure for generating 10 sample structures and evaluating the final value of AEyc as an 
average over these structures, the results of which are displayed in Table 6. For the neat ionic 
liquid systems, we see that solvation has little effect on the measure of AEyc as compared 
to AEhl for [pyrl4] [TFSI] and [pyrl3] [FSI ] . On the other hand, solvation systematically 
reduces the electrochemical window estimates for [EMIM]|BF 4 ], which could be related to 
the greater degree of aromaticity in the [EMIM] cation as compared to [pyrl4] and [pyrl3]. 
In all cases, the liquid phase computations result in the experimental electrochemical window 
being solidly bound by the pure and hybrid functionals, including [EMIM][BF 4 | in contrast 
to the cation/anion pair study. Additionally, the variation of AEyc between the neat and 
Li-doped species suggest that the addition of Li + has little influence on this measure of 
electrochemical window. For each system, the change in electrochemical window with varying 
degrees of Li-salt doping is on the order of 0.2 eV, which is approximately the uncertainty 
of our measures of AEyc (~ 0.2-0. 3). Essentially, this suggests the chemical effect of Li-salt 
doping is low as the cations and Li + have similar electrostatic influence on their surrounding 
solvent molecules. 


Conclusions 


The present work is an attempt to unite different levels of electronic and atomistic theory 
to provide an understanding of the physical properties of electrolytes for battery applications. 
We have combined zero-temperature DFT, DFT-MD, and PFF-MD to varying degrees to 
elucidate the structure, vibrations, transport, and electrochemical stability of three Li-doped 
ionic liquid electrolytes: [pyr 14] [TFSI] , [pyrl3][FSI], and [EMIM][BF 4 ], For [pyrl4] [TFSI] 
and [pyrl3][FSI], we find a Li + solvation structure that transitions from 2 -3 anions solvating 
Li + at low XLi to higher coordinated Li-network structures at Xl, > 0.1, which can conse- 
quently skew experimental measures of the number of anions in the Li + solvation shell. On 
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the other hand, [EMIM][BF 4 ] shows a tendency toward solvation shell with 4 coordinated 
anions at all levels of xl;. In general, fully solvated, charge neutral systems should be con- 
sidered to understand the favorability of a given Li + solvation number, However, it should 
be noted that cluster-level DFT computations successfully predict the most likely DFT-MD 
structures for a given value of n, and therefor represent a relatively inexpensive technique to 
survey likely solvation structures. Combining DFT-MD and PFF-MD at high-temperature, 
we have shown that DFT-MD is a viable technique for estimating the diffusion coefficient 
of Li + . PFF-MD suggests the majority of the error in the DFT-MD diffusion coefficients 
arises from the short time-scales used as opposed to size- related effects. Finally, we have 
determined that DFT measures of the lowest unoccupied and highest occupied electronic 
states may be a straightforward means of bounding the electrochemical window, with pure 
functionals being a lower bound and hybrid functionals being an upper bound. Solvation 
effects on the electrochemical window as evaluated through coupled DFT/PFF-MD simula- 
tion change the electrochemical window by only 5-10 %, except for the case of [EMIM]|BF 4 ], 
which exhibits a 20-25% change in electrochemical stability upon solvation. Furthermore, 
the effect of Li-salt doping on electrochemical window appears to be minor up to xl; = 0.33, 
with all lithiated measures being within the error of the neat value. 
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Table 1: Energetic decomposition of [Li + (Anion) n ]( n -!)- clusters as determined from 

B3LYP/6-31+G** calculations. Binding energy of the cluster (Eb), repulsive energy of 
the anions (Er), and attractive energy between the Li + and the anions (Ea) are given in 
kcal/mol. As a rule, the lowest energy clusters were produced when [TFSI] assumes the trans 
conformer and when [ESI] assumes the cis conformer in an rj 2 bond and trans in an z/ 1 bond. 




[Li(TFSI) n ] 

( n i) 

[Li(FSI) n ]( n 

-b- 

[Li(BF 4 ) n ] (n 

-b- 

n 

//-state 

Eb 

Er 

E a 

Eb 

Er 

e a 

Eb 

Er 

e a 

2 

2 z/ 2 

184.6 

60.1 

-244.7 

183.9 

61.0 

-244.9 

194.2 

75.7 

-269.9 

3 

CM 

CO 

146.5 

164.5 

-311.0 

146.2 

168.3 

-314.4 

147.8 

210.3 

-358.1 


2// 2 // 1 

149.2 

152.4 

-301.6 

148.1 

155.7 

-303.7 

151.3 

198.8 

-350.1 


'T]‘ 2 2ri 1 

150.8 

141.5 

-292.3 

148.9 

144.7 

-293.7 

t 




3// 1 

146.1 

133.9 

-280.0 

144.7 

135.7 

-280.3 

t 



4 

4-z/ 1 

81.6 

252.2 

-333.8 

79.4 

256.3 

-335.7 

59.5 

340.7 

-400.2 


f Configuration not stable 


34 



Table 2: Values of D app (f) (in units of 10~ 10 m 2 /s) for ions in ionic liquids at 363 K as 

evaluated using both DFT-MD and PFF-MD. The total MD simulation length ( t sim ) in 
picoseconds, the number of ion pairs in the MD simulation (N pa i rs ), and the time (t) used to 
evaluate D app (f) in picoseconds is given for each system. Values for the largest 144-216 N pairs 
PFF-MD systems are interpolated to the DFT-MD x t i values from independent simulations 
of xl; = 0.05, 0.10, and 0.15. Error as obtained from the standard deviation of the N° time 
origins used to average D app (f) is given for the DFT-MD results in parenthesis. 




^ sim 

N • 

1 ^ pairs 

t 

D+’ app (t) 

D-’ app (f) 

D Li ’ app (f) 

[pyrl4] [TFSI] 

DFT-MD 

0.1 

8 

0.05 

6.47 (3.2) 

6.10 (2.4) 

1.91 (0.6) 


PFF-MD 

6.0 

8 

0.05 

4.23 

3.17 

0.75 



6.0 

8 

t — > oo 

1.40 

0.89 

0.28 



6.0 

144 

t — > oo 

0.80 

0.66 

0.34 

[pyr 13] [FSI] 

DFT-MD 

0.1 

8 

0.05 

4.78 (2.3) 

4.27 (2.4) 

1.64 (0.5) 


PFF-MD 

6.0 

8 

0.05 

4.85 

3.85 

0.92 



6.0 

8 

t — > oo 

1.41 

1.29 

0.44 



6.0 

216 

t — > oo 

0.89 

0.91 

0.47 

[EMIM][BF 4 ] 

DFT-MD 

0.1 

8 

0.05 

7.77 (3.6) 

3.78 (2.2) 

1.03 (0.3) 


PFF-MD 

6.0 

8 

0.05 

4.77 

2.97 

0.97 



6.0 

8 

t — y oo 

1.97 

1.48 

0.75 



6.0 

216 

t — > oo 

1.93 

1.38 

0.65 
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Table 3: Room temperature measures of average Li + network size, (Mu.. .Li), the average 

number of anions in the solvation shell of a given Li + , (n), and total number of unique anions 
solvating all Li + , (N“), normalized by total number of Li-ions, Ny, as taken from PFF-MD 
simulation with 144-216 ion pairs. 



XLi 

(NLi...Li) 

( n ) 

(ND/Nl, 

[pyrl4] [TFSI] 

0.05 

1.1 

3.3 

3.3 


0.15 

1.4 

3.3 

2.8 


0.33 

1.5 

3.6 

2.2 

[pyrl3] [FSI] 

0.05 

1.1 

3.8 

3.7 


0.15 

1.4 

3.9 

3.1 


0.33 

1.6 

4.0 

2.5 

[EMIM][BF 4 ] 

0.05 

1.2 

3.9 

3.7 


0.15 

1.4 

3.9 

3.2 


0.33 

1.7 

4.0 

2.6 
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Table 4: Various DFT measures of the difference between the LUMO and HOMO energies 
(AEhl) and the sum of the ionization potential (IP) and electron affinity (EA) for diabatic 
(AEiE,d) and vertical (AEie, v ) excitations in single pairs of ionic liquids. The 6-31+G** 
basis set is used across all levels of theory and the units are eV. 




BPW91 

PW91 

PBE 

HSE06 

B3LYP 

M06 

[pyrl4] [TFSI] 

AE H l 

4.993 

4.781 

4.823 

6.267 

6.390 

6.340 


EA d 

0.052 

0.168 

0.134 

-0.017 

0.064 

-0.189 


IPd 

7.882 

8.141 

8.055 

8.442 

8.373 

8.875 


AEiE,d 

7.934 

8.309 

8.189 

8.425 

8.437 

8.686 


EA V 

-0.011 

0.117 

0.081 

-0.065 

0.014 

-0.247 


IPv 

8.719 

8.862 

8.779 

9.437 

9.356 

9.663 


AEie, v 

8.708 

8.979 

8.860 

9.372 

9.370 

9.426 

[pyrl3] [FSI] 

AE H l 

5.085 

4.823 

4.859 

6.364 

6.475 

6.479 


EA d 

2.150 

1.185 

1.078 

0.008 

0.099 

-0.183 


IPd 

8.178 

8.369 

8.288 

8.611 

8.569 

8.921 


AEiE,d 

10.328 

9.554 

9.366 

8.619 

8.668 

8.738 


EA V 

0.009 

0.177 

0.136 

-0.045 

0.039 

-0.232 


IPv 

9.109 

9.228 

9.148 

9.726 

9.660 

9.977 


AEie, v 

9.118 

9.405 

9.284 

9.681 

9.699 

9.745 

[EMIM][BF 4 ] 

AE H l 

4.549 

4.797 

4.660 

6.561 

6.841 

6.900 


EA d 

0.481 

0.436 

0.426 

0.338 

0.452 

0.243 


IPd 

8.284 

8.504 

8.423 

8.555 

8.608 

8.640 


AEiE,d 

8.765 

8.940 

8.849 

8.893 

9.060 

8.883 


EA V 

-0.096 

-0.027 

-0.058 

-0.261 

-0.164 

-0.386 


IPv 

9.551 

9.647 

9.576 

10.215 

10.214 

10.206 


AEie, v 

9.455 

9.620 

9.518 

9.954 

10.050 

9.820 
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Table 5: Influence of the number of ionic liquid cation/anions pairs (N pa i rs ) on liquid-phase 
DFT measures of the difference between the conduction band minimum (CBM) and valence 
band maximum (VBM) energies, referred to as AEyc- Energy is given in eV. 


A Eye 



N • 

PBE 

PW91 

HSE06 

B3LYP 

[pyrl4] [TFSI] 

12 

4.71 

4.71 

6.39 

6.66 


18 

4.58 

4.58 

6.25 

6.51 


24 

4.59 

4.59 

6.25 

6.53 

[pyrl3][FSI] 

12 

4.83 

4.86 

6.60 

6.86 


18 

4.71 

4.74 

6.46 

6.73 


24 

4.70 

4.73 

6.48 

6.74 

[EMIM][BF 4 ] 

12 

4.01 

4.00 

5.21 

5.96 


18 

3.91 

3.91 

5.15 

5.47 


24 

3.92 

3.92 

5.15 

5.46 
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Table 6: Influence of Li-salt mole fraction (xLi) on the difference in conduction band min- 
imum and valence band maximum energies (AEyc) as determined from liquid-phase DFT 
simulations having 24 ion pairs. Energies are given in units of eV. 


AEyc 



X Li + 

PBE 

HSE06 

B3LYP 

Exp. 

[pyrl4] [TFSI] 

0.000 

4.59 

6.25 

6.53 

3.8-6 36 ’ 37 ’ 66 


0.083 

4.61 

6.28 

6.57 



0.166 

4.50 

6.22 

6.51 



0.250 

4.50 

6.24 

6.53 



0.333 

4.36 

6.05 

6.34 


[pyrl3] [ESI] 

0.000 

4.70 

6.48 

6.74 

5.3 38 


0.083 

4.53 

6.32 

6.58 



0.166 

4.42 

6.25 

6.51 



0.250 

4.62 

6.43 

6.68 



0.333 

4.61 

6.44 

6.68 


[EMIM][BF 4 ] 

0.000 

3.92 

5.15 

5.46 

4 g 67,68 


0.083 

3.97 

5.22 

5.53 



0.166 

4.07 

5.30 

5.61 



0.250 

4.09 

5.31 

5.47 



0.333 

3.96 

5.18 

5.50 
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Figure 1: Energetically favorable [Li(Anion) n ]( n-1 )“ clusters from B3LYP/6-31+G** com- 
putations. Displayed are clusters having (a,d,g) n = 2, (b,e,h) n = 3, and (c,f,i) n = 4 anions 
for (a-c) [TFSI], (d-f) [FSI|, and (g-i) [BF 4 ], From DFT-MD simulations of liquids, solvation 
shells corresponding to clusters (a), (b), (d), (e) and (i) are found to be stable at 363 K. 
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Figure 2: Radial distribution function (g(r)) as obtained by PFF-MD and DFT-MD simu- 
lation at T — 363 K for Li + with the ionic liquid anions (a) [TFSI], (b) [FSI], and (c) [BF 4 ] 
as well as with the (d) O atoms in [TFSI], O atoms in [FSI], and F atoms in [BF4] . The 
radial distributions are averaged over a 100 ps DFT-MD trajectory and a 6 ns PFF-MD 
trajectory, with one Li + in ionic liquid systems having 8, 10, and 12 pairs for [pyrl4] [TFSI] , 
[pyrl3][FSI], and [EMIM][BF 4 ], respectively. 
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Figure 3: Raman active coupled CF 3 bend, 5 S (CF 3 ), SN stretch, z/ s (SN), modes of |TFSI] 
as determined by DFT for (a) n = 2, (b) n = 3, and (c) n = 4 [Li(TFSI) n ] ( n— x ) — clusters and 
(d) from the experiments of Lassegues and coworkers 9 for [BMIM]|TFSI| in both the neat 
form and that having xl; = 0.33 Li[TFSI|. 
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Figure 4: Raman active SN stretch, i/ s (SN), modes of [FSI] as determined by DFT for (a) n 
= 2, (b) n = 3, and (c) n = 4 [Li(FSI) n ]( n_1 )~ clusters and (d) from the experiments of Fujii 
and coworkers 12 for [EMIM][FSI| in both the neat form and that having = 0.225 Li[FSI] . 


43 



I I I I I I 

1210 1220 1230 1240 

to (cm’ 1 ) 


Figure 5: Raman active SO stretch, z/ s ( SO), modes of [FSI | as determined by DFT for (a) n 
= 2, (b) n = 3, and (c) n = 4 [Li(FSI) n ] < ' n ^ 1 0 clusters and (d) from the experiments of Fujii 
and coworkers 12 for [EMIM][FSI| in both the neat form and that having XLi = 0.225 Li[FSI|. 
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Figure 6: Infra-red signature of Li + in (a) [pyrl4] [TFSI] , (b) [pyrl3][FSI], and (c) 

[EMIM|[BF 4 ] as obtained from both 100 ps DFT-MD simulations at T = 363 K and DFT 
computations of the geometry optimized [Li(Anion) n ]( n "b~ clusters presented in Table 1. 
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Figure 7: Apparent diffusion coefficient of (D app (t)) as a function of time from (a-c) DFT- 
MD and (d-f) PFF-MD simulation at T — 363 K using small cells having 8, 10, and 12 pairs 
of [pyrl4][TFSI], [pyrl3][FSI], and [EMIM][BF 4 ], respectively, with one cation replaced by 
Li + . The values of D app as given are representative of a many time-origin averages, N° for 
both the 100 ps DFT-MD simulations (N° = 20) and the 6 ns PFF-MD simulations (N° = 
300). 
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